home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / a_utils / expanded.lha / expanded / src / Basis.C next >
Encoding:
C/C++ Source or Header  |  1992-03-19  |  25.9 KB  |  1,029 lines
  1. //
  2. // Linear-Affine-Projective Geometry Package
  3. //
  4. // Basis.C
  5. //
  6. // $Header$
  7. //
  8. // William J.R. Longabaugh 
  9. // University of Washington
  10. //
  11. // Implementation of the linear-affine-projective geometry
  12. // package described in William J.R. Longabaugh, "An Expanded
  13. // System for Coordinate-Free Geometric Programming", Master's 
  14. // thesis, University of Washington, 1992.
  15. //
  16. // Copyright (c) 1992, William J.R. Longabaugh
  17. //   Copying, use and development for non-commercial purposes permitted.
  18. //   All rights for commercial use reserved.
  19. //   This software is unsupported and without warranty; it is
  20. //   provided "as is".
  21. //
  22. // ***********************************************************************
  23.  
  24. #include <math.h>
  25. #include "Lap.h"
  26.  
  27. // ***********************************************************************
  28. // ***********************************************************************
  29. //
  30. // Basis Class
  31. //
  32. // ***********************************************************************
  33. // ***********************************************************************
  34. //
  35. // Private/protected member functions
  36. //
  37. // ***********************************************************************
  38. //
  39. // Used by other bases to build new bases:
  40. //
  41.  
  42. Basis::Basis(BasisType b, Space& ins,
  43.          char *namein, Matrix &m) : s(ins), tostdb(m), fromstdb(Inverse(m))
  44. {
  45.   type = ANY_BASIS;
  46.   holds = b;
  47.   strcpy(name, namein);
  48. }  
  49.  
  50. // ***********************************************************************
  51. //
  52. // Used to build a standard basis:
  53. //
  54.  
  55. Basis::Basis(BasisType b,
  56.          Space& ins, int size) : s(ins),
  57.                              tostdb(IdentityMatrix(size)),
  58.                      fromstdb(IdentityMatrix(size))
  59. {
  60.   type = ANY_BASIS;
  61.   holds = b;
  62.   strcpy(name, "Standard Basis");
  63. }  
  64.  
  65. // ***********************************************************************
  66. //
  67. // Dumps out Basis data:
  68. //
  69.  
  70. void Basis::data_out(ostream &c, int indent, char* n)
  71. {
  72.   char *ibloc = new char[indent + 1];
  73.   for (int i = 0; i < indent; i++) {
  74.     *(ibloc + i) = ' ';
  75.   }
  76.   *(ibloc + indent) = '\0';
  77.  
  78.   c << ibloc << ast;
  79.   c << ibloc << n;
  80.   c << ibloc << "Type of basis: ";
  81.   BasisTypeOut(c, type);
  82.   c << "\n";
  83.   c << ibloc << "Currently holds: ";
  84.   BasisTypeOut(c, holds);
  85.   c << "\n";
  86.   if (holds != NULL_BASIS) {
  87.     c << ibloc << "Name of basis: " << name << "\n";
  88.     c << ibloc << "Contained in space:\n";
  89.     s.debug_out(c, indent + ERR_IND);
  90.     c << ibloc << "Matrix representation of basis wrt standard basis:\n";
  91.     tostdb.debug_out(c, indent + ERR_IND);
  92.     c << ibloc << "Matrix representation of standard basis wrt this basis:\n";
  93.     fromstdb.debug_out(c, indent + ERR_IND);
  94.   }
  95.   c << ibloc << ast;
  96.  
  97.   delete ibloc;
  98.   return; 
  99. }
  100.  
  101. // ***********************************************************************
  102. //
  103. // Public member functions
  104. //
  105. // ***********************************************************************
  106. //
  107. // Need to do memberwise initialization, since Matrix class members need
  108. // to do memberwise initialization.
  109. //
  110.  
  111. Basis::Basis(Basis& b) : s(b.s), tostdb(b.tostdb), fromstdb(b.fromstdb)
  112. {
  113.   type = ANY_BASIS;
  114.   holds = b.holds;
  115.   b.Name(name);
  116. }  
  117.  
  118. // ***********************************************************************
  119. //
  120. // Need to do memberwise assignment, since Matrix class members need
  121. // to do memberwise assignment.
  122. //
  123.  
  124. Basis& Basis::operator=(Basis &b)
  125. {
  126. //
  127. // This unsavory checking is required to bypass the apparent inheritance of
  128. // assignment under g++ 1.37:
  129. //
  130.   if ((type != ANY_BASIS) &&
  131.       ((type != b.holds) && (b.holds != NULL_BASIS))) {
  132.     errh.ErrorExit("Basis& Basis::operator=(Basis&)",
  133.                    "Illegal assignment attempted", *this, b);
  134.   }
  135.   holds = b.holds;
  136.   s = b.s;
  137.   b.Name(name);
  138.   tostdb = b.tostdb;
  139.   fromstdb = b.fromstdb;
  140.   return (*this);
  141. }  
  142.  
  143. // ***********************************************************************
  144. //
  145. // Output for debugging:
  146. //
  147.  
  148. void Basis::debug_out(ostream &c, int indent)
  149. {
  150.   static char* name = "Basis object\n";
  151.   this->data_out(c, indent, name);
  152.   return; 
  153. }
  154.  
  155. // ***********************************************************************
  156. //
  157. // Return the nth member of the basis.
  158. //
  159.  
  160. GeOb Basis::operator[](int n)
  161. {
  162.   static char* sig = "GeOb Basis::operator[](int)";
  163.   int      maxn  = s.MatrixSize();
  164.   int       i;
  165.   GeObType rettype;
  166.   Space    retspace(s);
  167.   Matrix   hold;
  168.  
  169. // First check that the integer is in the correct range.  Note that
  170. // for projective spaces, the user can request the unit point of
  171. // the projective frame.
  172.  
  173.   if (holds == HFRAME) {
  174.     maxn++;
  175.   } 
  176.  
  177.   if ((n >= maxn) || (n < 0)) {
  178.     errh.ErrorExit(sig, "Specified n is invalid",
  179.                    ErrVal("Value of n = ", n), *this);
  180.   }
  181.  
  182. // The members of the basis are encoded in the matrix.  We need to
  183. // extract the correct row and build a GeOb using that row.  Two
  184. // special cases are the vector components of frames and the 
  185. // unit points of projective frames:
  186.  
  187.   if ((holds == FRAME) && (n != maxn - 1)) {
  188.     retspace = s.GetSpace(TANGENT);
  189.     rettype = AFF_VECTOR;
  190.     hold = (tostdb * Transpose(s.HoistTangent()))[n];
  191.   } else if ((holds == HFRAME) && (n == maxn - 1)) {
  192.     rettype = PROJ_POINT;
  193.     hold = tostdb[0];
  194.     for (i = 1; i < tostdb.Rows(); i++) {
  195.       hold += tostdb[i];
  196.     }
  197.   } else {  
  198.     rettype = s.NativeType();
  199.     hold = tostdb[n];
  200.   }
  201.  
  202.   GeOb retval(rettype, retspace, hold);
  203.   return (retval);
  204. }
  205.  
  206. // ***********************************************************************
  207. //
  208. // Apply the basis to the geometric object to get its coordinates
  209. // with respect to the basis.
  210. //
  211.  
  212. ScalarList Basis::operator()(GeOb &v)
  213. {
  214.   static char* sig = "ScalarList Basis::operator()(GeOb &)";
  215.   Boolean  cannot_map = FALSE;
  216.   GeOb       new_obj;
  217.   Matrix   retvals;
  218.   int       i;
  219.   GeObType native = s.NativeType();
  220.  
  221. // Check to see if object can be cast into the appropriate object
  222. // in the space of the basis, and do any necessary casting.  Then
  223. // get the coordinates with a matrix multiply.  Special case where
  224. // the basis is in an affine space and the object is from the
  225. // tangent space:
  226.  
  227.   if (v.CanMapTo(native)) {
  228.     new_obj = v.MapTo(native);
  229.     if (new_obj.SpaceOf() != s) {
  230.       cannot_map = TRUE;
  231.     } else {
  232.       retvals = new_obj.MatrixRep() * fromstdb;
  233.     }
  234.   } else if (((holds == SIMPLEX) || (holds == FRAME)) &&
  235.          v.CanMapTo(AFF_VECTOR)) {
  236.     new_obj = v.MapTo(AFF_VECTOR);
  237.     if (new_obj.SpaceOf() != s.GetSpace(TANGENT)) {
  238.       cannot_map = TRUE;
  239.     } else {
  240.       retvals = (new_obj.MatrixRep() * s.HoistTangent()) * fromstdb;
  241.     }
  242.   } else {
  243.     cannot_map = TRUE;
  244.   }
  245.  
  246. // Error exit if we cannot match spaces:
  247.  
  248.   if (cannot_map) {
  249.     errh.ErrorExit(sig, "Argument cannot be mapped to space of basis",
  250.                    v, *this);
  251.   }
  252.  
  253. // All that remains is to build a scalar list from the coordinates:
  254.  
  255.   ScalarList retlist(retvals.Columns());
  256.   for (i = 0; i < retvals.Columns(); i++) {
  257.     retlist[i] = retvals[0][i];
  258.   }
  259.   return (retlist);
  260. }
  261.  
  262. // ***********************************************************************
  263. //
  264. // Return the dual basis for a vector basis.
  265. //
  266.  
  267. VBasis Basis::Dual(void)
  268. {
  269.  
  270. // Must be a vector basis:
  271.  
  272.   if (holds != VBASIS) {
  273.     errh.ErrorExit("VBasis Basis::Dual(void)",
  274.                    "Basis is not a vector basis",
  275.                    *this);
  276.   }
  277.  
  278. // Copy the basis name, prepending with "Dual to".  If it is the
  279. // dual of a dual, return the original name:
  280.  
  281.   char buf[MAX_NAMESIZE];
  282.   char *dualto = "Dual to ";
  283.  
  284.   if (strncmp(name, dualto, 8) == 0) {
  285.     strcpy(buf, name + 8);
  286.     buf[MAX_NAMESIZE - 1] = '\0';
  287.   } else {
  288.     strncpy(buf, dualto, 8);
  289.     strncpy(buf + 8, name, MAX_NAMESIZE - 8);
  290.     buf[MAX_NAMESIZE - 1] = '\0';
  291.   }
  292.  
  293.   Basis retval(VBASIS, s.Dual(), buf, Transpose(Inverse(tostdb)));
  294.   return ((VBasis)retval);
  295. }
  296.  
  297. // ***********************************************************************
  298. // ***********************************************************************
  299. //
  300. // Simplex Class
  301. //
  302. // ***********************************************************************
  303. // ***********************************************************************
  304. //
  305. // Public member functions
  306. //
  307. // ***********************************************************************
  308. //
  309. // Need to do memberwise assignment, since Matrix class members need
  310. // to do memberwise assignment:
  311. //
  312.  
  313. Simplex& Simplex::operator=(Simplex &b)
  314. {
  315.   holds = b.holds;
  316.   s = b.s;
  317.   b.Name(name);
  318.   tostdb = b.tostdb;
  319.   fromstdb = b.fromstdb;
  320.   return (*this);
  321. }  
  322.  
  323. // ***********************************************************************
  324. //
  325. // Output for debugging:
  326. //
  327.  
  328. void Simplex::debug_out(ostream &c, int indent)
  329. {
  330.   static char* name = "Simplex object\n";
  331.   this->data_out(c, indent, name);
  332.   return; 
  333. }
  334.  
  335. // ***********************************************************************
  336. //
  337. // Works if the GeOb is a tuple of points from a cartesian product
  338. // affine space (A x A x A x ...)
  339. //
  340.  
  341. Simplex::Simplex(char* namein, GeOb &t)
  342. {
  343.   static char* sig = "Simplex::Simplex(char*, GeOb&)";
  344.   int i;
  345.   int tuplesize;
  346.   GeOb new_obj;
  347.   SpaceType tupletype;
  348.  
  349. // Set the type:
  350.  
  351.   type = SIMPLEX;
  352.   holds = SIMPLEX;
  353.  
  354. // Local copy of the debug name:
  355.  
  356.   strncpy(name, namein, MAX_NAMESIZE - 1);
  357.   name[MAX_NAMESIZE - 1] = '\0';
  358.  
  359. // Check to see if the space is a CP Affine space:
  360.  
  361.   tuplesize = (t.SpaceOf()).CPSpaceSize();
  362.   tupletype = (t.SpaceOf()).Holds();
  363.  
  364.   if ((tupletype != AFF_SPACE) || (tuplesize == 1)) {
  365.     errh.ErrorExit(sig, "Space is not a cartesian product affine space", t);
  366.   }
  367.  
  368.   s = t.SpaceOf()[0];
  369.   tostdb = Matrix(s.MatrixSize());
  370.  
  371. // Since the point tuple is from an affine space, each tuple element
  372. // is an affine point.  We just need to check that the tuple is the
  373. // correct size and that all points are from the same space as we build
  374. // the matrix:
  375.  
  376.   if (tuplesize != s.MatrixSize()) {
  377.     errh.ErrorExit(sig, "Tuple is not the correct size for the space", t, s);
  378.   }
  379.  
  380.   for (i = 0; i < s.MatrixSize(); i++) {
  381.     new_obj = t[i];
  382.     if (new_obj.SpaceOf() != s) {
  383.       errh.ErrorExit(sig, "Tuple element not from correct space", new_obj, s);
  384.     }
  385.     tostdb[i] = (new_obj.MatrixRep())[0];
  386.   }
  387.  
  388. // Make sure that the simplex is valid by checking that all the
  389. // points are independent.  If it is OK, invert it:
  390.  
  391.   if (fabs(Det(tostdb)) < EPSILON) {
  392.       errh.ErrorExit(sig, "Tuple elements are not independent", t);
  393.   }
  394.   fromstdb = Inverse(tostdb);
  395. }
  396.  
  397. // ***********************************************************************
  398. //
  399. // Builds a Simplex in the specified space:
  400. //
  401.  
  402. Simplex::Simplex(char* namein, ASpace& ins, GeObList& t)
  403. {
  404.   static char* sig = "Simplex::Simplex(char*, Space&, GeObList&)";
  405.   int i;
  406.   Boolean cannot_map = FALSE;
  407.   GeOb new_obj;
  408.  
  409. // Set the type:
  410.  
  411.   type = SIMPLEX;
  412.   holds = SIMPLEX;
  413.  
  414. // Local copy of the debug name:
  415.  
  416.   strncpy(name, namein, MAX_NAMESIZE - 1);
  417.   name[MAX_NAMESIZE - 1] = '\0';
  418.  
  419.   s = ins;
  420.   tostdb = Matrix(s.MatrixSize());
  421.  
  422. // Check to see if we have the correct number of points, and build
  423. // the basis matrix using the points:
  424.  
  425.   if (t.Length() != s.MatrixSize()) {
  426.     errh.ErrorExit(sig, "Incorrect number of objects specified", t, ins);
  427.   }
  428.  
  429.   for (i = 0; i < s.MatrixSize(); i++) {
  430.     cannot_map = FALSE;
  431.     if (t[i].CanMapTo(AFF_POINT)) {
  432.         new_obj = t[i].MapTo(AFF_POINT);
  433.         if (new_obj.SpaceOf() != s) {
  434.       cannot_map = TRUE;
  435.     } else {
  436.        tostdb[i] = (new_obj.MatrixRep())[0];
  437.     }
  438.     } else {
  439.       cannot_map = TRUE;
  440.     }
  441.     if (cannot_map) {
  442.       errh.ErrorExit(sig, "Object cannot be mapped to specified space",
  443.                      t[i], s);
  444.     }
  445.   }
  446.  
  447. // Make sure that the simplex is valid by checking that all the
  448. // points are independent.  If it is OK, invert it:
  449.  
  450.   if (fabs(Det(tostdb)) < EPSILON) {
  451.       errh.ErrorExit(sig, "Points are not independent", t);
  452.   }
  453.   fromstdb = Inverse(tostdb);
  454. }
  455.  
  456. // ***********************************************************************
  457. //
  458. // Used to cast a general basis down to a Simplex.  This is also the
  459. // way to change a Frame into a Simplex
  460. //
  461.  
  462. Simplex::Simplex(Basis &f) : (f)
  463. {
  464.   static char* sig = "Simplex::Simplex(Basis&)";
  465.   int i;
  466.  
  467.   type = SIMPLEX;
  468.  
  469. // Conversion of Frames to Simplices:
  470.  
  471.   if (holds == FRAME) {
  472.  
  473.     holds = SIMPLEX;
  474.  
  475. //  Add a prefix to the name:
  476.  
  477.     char buf[MAX_NAMESIZE + 13];
  478.     char *smplxfrom = "Simplex from ";
  479.     strncpy(buf, smplxfrom, 13);
  480.     this->Name(buf + 13);
  481.     strncpy(name, buf, MAX_NAMESIZE - 1);
  482.     name[MAX_NAMESIZE - 1] = '\0';
  483.  
  484. //  Need to build the new matrix by taking the frame origin and
  485. //  adding it to all the vectors to make points.  The origin stays
  486. //  as the last row of the target matrix.
  487.  
  488.     for (i = 0; i < (s.MatrixSize() - 1); i++) {
  489.       tostdb[i] = (tostdb[s.MatrixSize() - 1] + tostdb[i])[0];
  490.     }
  491.  
  492. //  Make sure that the simplex is valid.  If it is OK, invert it:
  493.  
  494.     if (fabs(Det(tostdb)) < EPSILON) {
  495.       errh.ErrorExit(sig, "Points are not independent", *this);
  496.     }
  497.     fromstdb = Inverse(tostdb);
  498.  
  499.   } else if (holds != SIMPLEX) {
  500.     errh.ErrorExit(sig, "Attempted to cast a non-Simplex to a Simplex", f);
  501.   }
  502. }
  503.  
  504. // ***********************************************************************
  505. // ***********************************************************************
  506. //
  507. // Frame Class
  508. //
  509. // ***********************************************************************
  510. // ***********************************************************************
  511. //
  512. // Public member functions
  513. //
  514. // ***********************************************************************
  515. //
  516. // Used to cast a general basis down to a Frame.
  517. //
  518.  
  519. Frame::Frame(Basis &f) : (f)
  520. {
  521.   type = FRAME;
  522.   if ((holds != FRAME) && (holds != NULL_BASIS)) {
  523.       errh.ErrorExit("Frame::Frame(Basis &)",
  524.                      "Attempted to cast a non-frame to a frame", f);
  525.   }
  526. }
  527.  
  528. // ***********************************************************************
  529. //
  530. // Need to do memberwise assignment, since Matrix class members need
  531. // to do memberwise assignment:
  532. //
  533.  
  534. Frame& Frame::operator=(Frame &f)
  535. {
  536.   holds = f.holds;
  537.   s = f.s;
  538.   f.Name(name);
  539.   tostdb = f.tostdb;
  540.   fromstdb = f.fromstdb;
  541.   return (*this);
  542. }  
  543.  
  544. // ***********************************************************************
  545. //
  546. // Output for debugging:
  547. //
  548.  
  549. void Frame::debug_out(ostream &c, int indent)
  550. {
  551.   static char* name = "Frame object\n";
  552.   this->data_out(c, indent, name);
  553.   return; 
  554. }
  555.  
  556. // ***********************************************************************
  557. //
  558. // Builds a Frame in the specified space:
  559. //
  560.  
  561. Frame::Frame(char* namein, ASpace& ins, GeObList& t)
  562. {
  563.   static char* sig = "Frame::Frame(char*, Space&, GeObList&)";
  564.   int i;
  565.   Boolean cannot_map = FALSE;
  566.   GeOb new_obj;
  567.  
  568. // Set the type:
  569.  
  570.   type = FRAME;
  571.   holds = FRAME;
  572.  
  573. // Local copy of the debug name:
  574.  
  575.   strncpy(name, namein, MAX_NAMESIZE - 1);
  576.   name[MAX_NAMESIZE - 1] = '\0';
  577.  
  578.   s = ins;
  579.   tostdb = Matrix(s.MatrixSize());
  580.  
  581. // Check to see if we have the correct number of objects, and build
  582. // the basis matrix using the objects:
  583.  
  584.   if (t.Length() != s.MatrixSize()) {
  585.     errh.ErrorExit(sig, "Incorrect number of objects specified", t, ins);
  586.   }
  587.  
  588. // Map all but the last object to vectors:
  589.  
  590.   for (i = 0; (i < (s.MatrixSize() - 1)) && (cannot_map == FALSE); i++) {
  591.     if (t[i].CanMapTo(AFF_VECTOR)) {
  592.         new_obj = t[i].MapTo(AFF_VECTOR);
  593.         if (new_obj.SpaceOf() != s.GetSpace(TANGENT)) {
  594.       cannot_map = TRUE;
  595.     } else {
  596.        tostdb[i] = (new_obj.MatrixRep() * s.HoistTangent())[0];
  597.     }
  598.     } else {
  599.       cannot_map = TRUE;
  600.     }
  601.   }
  602.  
  603. // Map last object to a point:
  604.  
  605.   if (t[s.MatrixSize() - 1].CanMapTo(AFF_POINT)) {
  606.     new_obj = t[s.MatrixSize() - 1].MapTo(AFF_POINT);
  607.     if (new_obj.SpaceOf() != s) {
  608.       cannot_map = TRUE;
  609.     } else {
  610.       tostdb[s.MatrixSize() - 1] = (new_obj.MatrixRep())[0];
  611.     }
  612.   } else {
  613.     cannot_map = TRUE;
  614.   }
  615.  
  616.   if (cannot_map) {
  617.     errh.ErrorExit(sig, "Object cannot be mapped to specified space", t, s);
  618.   }
  619.  
  620. // Make sure that the frame is valid by checking that all the
  621. // vectors are independent.  If it is OK, invert it:
  622.  
  623.   if (fabs(Det(tostdb)) < EPSILON) {
  624.       errh.ErrorExit(sig, "Vectors are not independent", t);
  625.   }
  626.   fromstdb = Inverse(tostdb);
  627. }
  628.  
  629. // ***********************************************************************
  630. //
  631. // Used to convert a Simplex into a Frame by selecting a
  632. // point to serve as the origin:
  633. //
  634.  
  635. Frame::Frame(Simplex& b, int n)
  636. {
  637.   static char* sig = "Frame::Frame(Simplex&, int)";
  638.   int i;
  639.   int slot;
  640.  
  641. // Set the type, name, and space:
  642.  
  643.   type = FRAME;
  644.   holds = FRAME;
  645.   s = b.SpaceOf();
  646.   tostdb = Matrix(s.MatrixSize());
  647.  
  648. // Prefix the name before copying:
  649.  
  650.   char buf[MAX_NAMESIZE + 11];
  651.   char *framefrom = "Frame from ";
  652.   strncpy(buf, framefrom, 11);
  653.   b.Name(buf + 11);
  654.   strncpy(name, buf, MAX_NAMESIZE - 1);
  655.   name[MAX_NAMESIZE - 1] = '\0';
  656.  
  657. // Check to see if the specified n is legal:
  658.  
  659.   if ((n < 0) || (n >= s.MatrixSize())) {
  660.     errh.ErrorExit(sig, "Specified n is out of range",
  661.            ErrVal("n = ", n), b);
  662.   }
  663.  
  664. // Need to build the new matrix by taking the chosen point and
  665. // subtracting it from all the the other points to make vectors.
  666. // The chosen point is put in the last row of the target matrix.
  667.  
  668.  
  669. // Map all but the last object to vectors:
  670.  
  671.   for (slot = 0, i = 0; i < s.MatrixSize(); i++) {
  672.     if (i == n) {
  673.       tostdb[s.MatrixSize() - 1] = (b.Tostdb())[n];
  674.     } else {
  675.       tostdb[slot++] = ((b.Tostdb())[i] - (b.Tostdb())[n])[0];
  676.     }
  677.   }
  678.  
  679. // Make sure that the frame is valid.  If it is OK, invert it:
  680.  
  681.   if (fabs(Det(tostdb)) < EPSILON) {
  682.       errh.ErrorExit(sig, "Vectors are not independent", *this);
  683.   }
  684.   fromstdb = Inverse(tostdb);
  685. }
  686.  
  687. // ***********************************************************************
  688. // ***********************************************************************
  689. //
  690. // VBasis Class
  691. //
  692. // ***********************************************************************
  693. // ***********************************************************************
  694. //
  695. // Public member functions
  696. //
  697. // ***********************************************************************
  698. //
  699. // Used to cast a general basis down to a VBasis.
  700. //
  701.  
  702. VBasis::VBasis(Basis& f) : (f)
  703. {
  704.   type = VBASIS;
  705.   if ((holds != VBASIS) && (holds != NULL_BASIS)) {
  706.       errh.ErrorExit("VBasis::VBasis(Basis&)",
  707.                      "Attempted to cast a non-VBasis to a VBasis", f);
  708.   }
  709. }
  710.  
  711. // ***********************************************************************
  712. //
  713. // Need to do memberwise assignment, since Matrix class members need
  714. // to do memberwise assignment:
  715. //
  716.  
  717. VBasis& VBasis::operator=(VBasis& b)
  718. {
  719.   holds = b.holds;
  720.   s = b.s;
  721.   b.Name(name);
  722.   tostdb = b.tostdb;
  723.   fromstdb = b.fromstdb;
  724.   return (*this);
  725. }  
  726.  
  727. // ***********************************************************************
  728. //
  729. // Output for debugging:
  730. //
  731.  
  732. void VBasis::debug_out(ostream &c, int indent)
  733. {
  734.   static char* name = "VBasis Object\n";
  735.   this->data_out(c, indent, name);
  736.   return; 
  737. }
  738.  
  739. // ***********************************************************************
  740. //
  741. // Works if the GeOb is a tuple of vectors from a cartesian product
  742. // vector space.
  743. //
  744.  
  745. VBasis::VBasis(char* namein, GeOb& t)
  746. {
  747.   static char* sig = "VBasis::VBasis(char*, GeOb&)";
  748.   int       i;
  749.   int       tuplesize;
  750.   GeOb      new_obj;
  751.   SpaceType tupletype;
  752.   Space     firstspace;
  753.   GeObType  targettype;
  754.   Boolean   cannot_map = FALSE;
  755.   
  756.  
  757. // Set the type:
  758.  
  759.   type = VBASIS;
  760.   holds = VBASIS;
  761.  
  762. // Local copy of the debug name:
  763.  
  764.   strncpy(name, namein, MAX_NAMESIZE - 1);
  765.   name[MAX_NAMESIZE - 1] = '\0';
  766.  
  767. // Check to see if the space is a CP Vector space:
  768.  
  769.   tuplesize = (t.SpaceOf()).CPSpaceSize();
  770.   tupletype = (t.SpaceOf()).Holds();
  771.  
  772.   if ((tupletype != VEC_SPACE) || (tuplesize == 1)) {
  773.     errh.ErrorExit(sig, "Space is not a cartesian product vector space", t);
  774.   }
  775.  
  776. // Figure out the space to build the basis in.  The vector tuple could
  777. // be a mixture of vectors from a linearization/tangent space pair.
  778. // Which space to use depends on the size of the tuple and the dimensions
  779. // of the spaces:
  780.  
  781.   firstspace = t.SpaceOf()[0];
  782.   if (firstspace.Dim() == tuplesize) {
  783.     s = firstspace;
  784.   } else if (firstspace.Dim() == tuplesize + 1) {
  785.     s = firstspace.GetSpace(TANGENT);
  786.   } else if (firstspace.Dim() == tuplesize - 1) {
  787.     s = firstspace.GetSpace(LINEARIZATION);
  788.   } else {
  789.     errh.ErrorExit(sig, "Tuple is not the correct size for the space",
  790.            t, firstspace);
  791.   }
  792.    
  793.   targettype = s.NativeType();
  794.   tostdb = Matrix(s.MatrixSize());
  795.  
  796. // Now build up the matrix by mapping the vectors into the desired space: 
  797.  
  798.   for (i = 0; i < s.MatrixSize(); i++) {
  799.     cannot_map = FALSE;
  800.     if (t[i].CanMapTo(targettype)) {
  801.         new_obj = t[i].MapTo(targettype);
  802.         if (new_obj.SpaceOf() != s) {
  803.       cannot_map = TRUE;
  804.     } else {
  805.        tostdb[i] = (new_obj.MatrixRep())[0];
  806.     }
  807.     } else {
  808.       cannot_map = TRUE;
  809.     }
  810.     if (cannot_map) {
  811.       errh.ErrorExit(sig, "Tuple element cannot be mapped to specified space",
  812.                      t[i], s);
  813.     }
  814.   }
  815.  
  816. // Make sure that the basis is valid by checking that all the
  817. // vectors are independent.  If it is OK, invert it:
  818.  
  819.   if (fabs(Det(tostdb)) < EPSILON) {
  820.       errh.ErrorExit(sig, "Vectors are not independent", t);
  821.   }
  822.   fromstdb = Inverse(tostdb);
  823. }
  824.  
  825. // ***********************************************************************
  826. //
  827. // Builds a VBasis in the specified space:
  828. //
  829.  
  830. VBasis::VBasis(char* namein, VSpace& ins, GeObList& t)
  831. {
  832.   static char* sig = "VBasis::VBasis(char*, Space&, GeObList&)";
  833.   int i;
  834.   Boolean cannot_map = FALSE;
  835.   GeOb new_obj;
  836.   GeObType targettype;
  837.  
  838. // Set the type:
  839.  
  840.   type = VBASIS;
  841.   holds = VBASIS;
  842.  
  843. // Local copy of the debug name:
  844.  
  845.   strncpy(name, namein, MAX_NAMESIZE - 1);
  846.   name[MAX_NAMESIZE - 1] = '\0';
  847.  
  848.   s = ins;
  849.   tostdb = Matrix(s.MatrixSize());
  850.  
  851. // Check to see if we have the correct number of vectors, and build
  852. // the basis matrix using the vectors:
  853.  
  854.   if (t.Length() != s.MatrixSize()) {
  855.     errh.ErrorExit(sig, "Incorrect number of objects specified", t, ins);
  856.   }
  857.  
  858. // Map list elements to correct types:
  859.  
  860.   targettype = s.NativeType();
  861.  
  862.   for (i = 0; i < s.MatrixSize(); i++) {
  863.     cannot_map = FALSE;
  864.     if (t[i].CanMapTo(targettype)) {
  865.         new_obj = t[i].MapTo(targettype);
  866.         if (new_obj.SpaceOf() != s) {
  867.       cannot_map = TRUE;
  868.     } else {
  869.        tostdb[i] = (new_obj.MatrixRep())[0];
  870.     }
  871.     } else {
  872.       cannot_map = TRUE;
  873.     }
  874.     if (cannot_map) {
  875.       errh.ErrorExit(sig, "Object cannot be mapped to specified space",
  876.                      t[i], s);
  877.     }
  878.   }
  879.  
  880. // Make sure that the basis is valid by checking that all the
  881. // vectors are independent.  If it is OK, invert it:
  882.  
  883.   if (fabs(Det(tostdb)) < EPSILON) {
  884.       errh.ErrorExit(sig, "Vectors are not independent", t);
  885.   }
  886.   fromstdb = Inverse(tostdb);
  887. }
  888.  
  889. // ***********************************************************************
  890. // ***********************************************************************
  891. //
  892. // HFrame Class
  893. //
  894. // ***********************************************************************
  895. // ***********************************************************************
  896. //
  897. // Public member functions
  898. //
  899. // ***********************************************************************
  900. //
  901. // Used to cast a general basis down to an HFrame.
  902. //
  903.  
  904. HFrame::HFrame(Basis& f) : (f)
  905. {
  906.   type = HFRAME;
  907.   if ((holds != HFRAME) && (holds != NULL_BASIS)) {
  908.       errh.ErrorExit("HFrame::HFrame(Basis &)",
  909.                      "Attempted to cast a non-HFrame to an HFrame", f);
  910.   }
  911. }
  912.  
  913. // ***********************************************************************
  914. //
  915. // Need to do memberwise assignment, since Matrix class members need
  916. // to do memberwise assignment:
  917. //
  918.  
  919. HFrame& HFrame::operator=(HFrame& b)
  920. {
  921.   holds = b.holds;
  922.   s = b.s;
  923.   b.Name(name);
  924.   tostdb = b.tostdb;
  925.   fromstdb = b.fromstdb;
  926.   return (*this);
  927. }  
  928.  
  929. // ***********************************************************************
  930. //
  931. // Output for debugging:
  932. //
  933.  
  934. void HFrame::debug_out(ostream &c, int indent)
  935. {
  936.   static char* name = "HFrame Object\n";
  937.   this->data_out(c, indent, name);
  938.   return; 
  939. }
  940.  
  941. // ***********************************************************************
  942. //
  943. // Builds a Projective Frame in the specified space:
  944. //
  945.  
  946. HFrame::HFrame(char* namein, PSpace& ins, GeObList& t)
  947. {
  948.   static char* sig = "HFrame::HFrame(char*, Space&, GeObList&)";
  949.   int i;
  950.   Boolean cannot_map = FALSE;
  951.   GeOb new_obj;
  952.   Matrix unit_pt, tempinv, coeff;
  953.  
  954. // Set the type:
  955.  
  956.   type = HFRAME;
  957.   holds = HFRAME;
  958.  
  959. // Local copy of the debug name:
  960.  
  961.   strncpy(name, namein, MAX_NAMESIZE - 1);
  962.   name[MAX_NAMESIZE - 1] = '\0';
  963.  
  964.   s = ins;
  965.   tostdb = Matrix(s.MatrixSize());
  966.  
  967. // Check to see if we have the correct number of points, and build
  968. // up a temporary matrix and a unit point matrix using the points.  
  969.  
  970.   if (t.Length() != (s.MatrixSize() + 1)) {
  971.     errh.ErrorExit(sig, "Incorrect number of objects specified", t, ins);
  972.   }
  973.  
  974.   Matrix tempmtx(s.MatrixSize());
  975.  
  976.   for (i = 0; (i <= s.MatrixSize()) && (cannot_map == FALSE); i++) {
  977.     cannot_map = FALSE;
  978.     if (t[i].CanMapTo(PROJ_POINT)) {
  979.         new_obj = t[i].MapTo(PROJ_POINT);
  980.         if (new_obj.SpaceOf() != s) {
  981.       cannot_map = TRUE;
  982.     } else if (i == s.MatrixSize()) {
  983.        unit_pt = new_obj.MatrixRep();
  984.     } else {
  985.        tempmtx[i] = (new_obj.MatrixRep())[0];
  986.     }
  987.     } else {
  988.       cannot_map = TRUE;
  989.     }
  990.   }
  991.  
  992.   if (cannot_map) {
  993.     errh.ErrorExit(sig, "Object cannot be mapped to specified space",
  994.                    t[i], s);
  995.   }
  996.  
  997. // We need to scale the rows of the temporary matrix so that the
  998. // unit point has homogenous coordinates (1, 1, 1, ...).  Do this
  999. // by solving a system of equations.  Invert the temporary matrix:
  1000.  
  1001.   if (fabs(Det(tempmtx)) < EPSILON) {
  1002.     errh.ErrorExit(sig, "Points are not independent", t);
  1003.   }
  1004.   tempinv = Inverse(tempmtx);
  1005.   
  1006. // Solve for the row coefficients.  If any coefficient = 0, the
  1007. // unit point was not in general position.  Otherwise, multiply
  1008. // the matrix rows by the coefficients:
  1009.  
  1010.   coeff = unit_pt * tempinv;
  1011.   for (i = 0; i < s.MatrixSize(); i++) {
  1012.     if (fabs(coeff[0][i]) < EPSILON) {
  1013.       errh.ErrorExit(sig, "Unit point not in general position", t);
  1014.     } else {
  1015.       tostdb[i] = (coeff[0][i] * tempmtx[i])[0];
  1016.     }
  1017.   }
  1018.  
  1019. // Make sure that the hframe is valid by checking again that the
  1020. // matrix is invertible.  If it is OK, invert it:
  1021.  
  1022.   if (fabs(Det(tostdb)) < EPSILON) {
  1023.     errh.ErrorExit(sig, "Points are not independent", t);
  1024.   }
  1025.   fromstdb = Inverse(tostdb);
  1026. }
  1027.  
  1028. // ***********************************************************************
  1029.